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Abstract — The  identification  of  a  dynamic,  nonlinear  model  of  human 
ankle  stiffness  is  considered  in  a  minimum  mean  squared  error  frame¬ 
work.  The  model  consists  of  two  parallel  pathways,  one  representing  the 
intrinsic  dynamics,  the  other  representing  the  reflex  contribution  to  the 
stiffness.  The  model  is  shown  to  be  linear  in  all  of  its  parameters,  except 
for  those  used  to  describe  a  single  static  nonlinearity  in  the  reflex  path¬ 
way.  A  separable  least  squares  optimization  algorithm  is  developed  which 
takes  advantage  of  this  structure.  This  new  algorithm  is  applied  to  ex¬ 
perimental  stretch  reflex  data,  and  the  results  compared  to  the  current 
state-of-the-art  algorithm,  an  iterative  technique  which  fits  the  two  path¬ 
ways  alternately.  The  relative  merits  of  the  two  approaches  are  discussed. 

Keywords —  nonlinear  system  identification,  mean  squared  optimiza¬ 
tion,  separable  least  squares,  stretch  reflex  dynamics. 


I.  Introduction 

System  identification  methods,  which  construct  mathemat¬ 
ical  models  of  dynamic  systems  from  measurements  of  their 
inputs  and  outputs,  are  often  used  to  study  physiological  sys¬ 
tems.  Since  these  systems  are  often  highly  nonlinear,  physio¬ 
logical  applications  often  require  nonlinear  system  identifica¬ 
tion  methods  [9]. 

One  of  the  major  advances  in  the  field  of  nonlinear  system 
identification  has  been  the  adoption  of  explicit  least  squares 
estimation  methods  in  place  of  earlier  cross-correlation  based 
algorithms  [9].  For  models,  like  the  Volterra  and  Wiener  se¬ 
ries,  which  are  linear  in  their  parameters,  the  minimum  mean 
squared  error  (MMSE)  solution  can  be  found  in  closed  form, 
simply  by  solving  a  linear  regression.  However,  the  number 
of  parameters  required  to  represent  high-order  kernels  limits 
these  methods  to  relatively  low-order  systems  (second  or  per¬ 
haps  third  order  nonlinearities). 

Systems  that  include  high-order  nonlinearities  are  often  rep¬ 
resented  using  block-structured  models:  interconnections  of 
zero-memory  (i.e.  static)  nonlinearities  and  dynamic  linear 
systems.  The  simplest  of  these  structures  are  the  Wiener  mod¬ 
el,  a  linear  dynamic  element  followed  by  a  static  nonlinearity, 
and  the  Hammerstein  model,  a  static  nonlinearity  followed  by 
a  linear  filter  [2], 

Block  structured  models  are  not  very  general.  Thus,  the 
topology  of  the  model  must  be  appropriate  for  the  system  being 
studied.  Figure  1  shows  a  block  structured  model  that  is  used 
to  represent  the  dynamic  stiffness  of  the  human  ankle  [5].  The 
upper  pathway  represents  the  “intrinsic  stiffness”  of  the  ankle, 
the  dynamic  relationship  between  its  position  and  the  result¬ 
ing  torque,  in  the  absence  of  any  reflexive  actions.  The  lower 
pathway  represents  the  contribution  of  the  stretch  reflex  to  the 
overall  stiffness.  The  position  input  is  differentiated,  produc¬ 


ing  the  angular  velocity  of  the  ankle.  This  velocity  signal  is 
processed  by  a  static  nonlinearity,  which  may  represent  neural 
encoding  processes  in  the  muscle  spindle  and  in  the  alpha  mo¬ 
toneuron  pool.  The  final  element  in  this  pathway  is  a  dynamic 
linear  system,  which  is  thought  to  represent  the  contractile  dy¬ 
namics  of  the  muscle. 


Intrinsic  Stiffness  Pathway 


Reflex  Stiffness  Pathway 


Fig.  1.  Parallel  Cascade  model  of  the  reflex  stiffness,  the  dynamic  relationship 
between  the  ankle  angle,  6,  and  the  ankle  torque,  Tq.  The  model  includes 
both  intrinsic  (upper  pathway)  and  reflex  (lower  pathway)  components. 


Like  the  other  block-structured  models,  this  parallel  cas¬ 
cade  stiffness  model  (PCSM)  isn’t  linear  in  its  parameters. 
Thus,  unlike  the  Volterra  series,  there  is  no  generally  applica¬ 
ble  closed-form  solution  for  the  optimal  (MMSE)  parameters. 
However,  like  the  Hammerstein  system  [10],  the  PC  stiffness 
model  is  only  nonlinear  in  the  parameters  describing  its  stat¬ 
ic  nonlinearity,  and  linear  in  the  remaining  parameters.  Thus, 
this  structure  is  an  ideal  candidate  for  an  identification  tech¬ 
nique  based  on  separable  least  squares  optimization. 

In  this  paper,  we  will  develop  an  identification  technique 
for  the  PCSM  based  on  a  separable  least  squares  optimization. 
This  new  technique  will  be  compared  to  the  current  state  of  the 
art,  an  iterative,  correlation  based  algorithm,  using  experimen¬ 
tal  stretch  reflex  data.  The  physiological  insights  gained  from 
the  newly  identified  models  will  then  be  discussed. 

II.  Theory 

The  input  and  output  of  the  PCSM,  shown  in  Fig.  1,  are  the 
ankle  (angular)  position,  9(t ),  and  the  torque,  Tq(t),  generated 
about  it.  The  upper  pathway  represents  the  intrinsic  stiffness. 
Its  output  is  the  torque  that  would  be  produced  in  the  absence 
of  reflexes.  The  transfer  function  representing  dynamic  stiff¬ 
ness  is  improper,  in  that  it  has  more  zeros  than  poles  [4],  In 
the  discrete-time  domain,  this  means  that  the  impulse  response 
of  the  intrinsic  stiffness  will  be  a  two-sided  filter  [1],  Thus, 
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the  intrinsic  stiffness  is  represented  by  the  two-sided  impulse 
response  (IRF). 

Li 

Ti(t)=  ^2  hiS(T)0(t  —  t)  (1) 

r-i-  /./ 

The  lower  pathway  represents  the  reflex  contribution  to  the 
total  ankle  stiffness.  The  stretch  reflex  EMG  is  often  modeled 
as  a  Hammerstein  system  [10]  with  ankle  velocity  as  its  input. 
In  the  PCSM,  the  lower  pathway  consists  of  a  differentiator, 
which  computes  the  ankle  velocity,  followed  by  a  Hammer¬ 
stein  system,  which  represents  the  stretch  reflex. 

Let  9(t)  be  the  angular  velocity,  obtained  by  numerically 
differentiating  the  position,  6(t ),  and  let  the  output  of  the  static 
nonlinearity,  x(t),  be  computed  using  a  degree  Q  polynomial, 

Q 

x(t)  =  '^2cqdq(t)  (2) 

9=0 

The  linear  subsystem  in  the  reflex  pathway  will  be  represented 
by  its  impulse  response. 

Lr 

TR{t)  =  hRs(r)x(t  -  t)  (3) 

T— 0 

which,  unlike  the  IRF  of  the  intrinsic  stiffness,  is  causal.  The 
PCSM  output  is  obtained  by  summing  Tp>s(t)  and  Tis(t), 

Li 

Tq(t)  =  ^2  his(r)9(t-T)  (4) 

r=-Li 

Lr  (  Q  \ 

+  hRs{r)  ~  T) 

r=0  \q—0  ) 

A.  Iterative,  Correlation-Based  Identification 

Currently,  the  elements  of  the  PCSM  are  identified  using 
an  iterative,  cross-correlation  based  scheme.  This  is  possible 
because  of  the  delays  inherent  in  the  system.  The  two-sided 
impulse  response  representing  the  inherent  dynamics  dies  out 
well  before  40  ms.  The  reflex  component,  on  the  other  hand, 
includes  a  delay  of  at  least  40  ms,  due  to  propagation  delays. 
Thus,  the  two  contributions  can  be  separated  temporally. 

The  first  step  in  fitting  the  reflex  cascade  model  is  to  fit 
a  two-sided  IRF,  /i/s(t),  whose  causal  dynamics  are  limited 
to  ±40?ns,  between  9(t),  and  Tq(t)  using  a  cross-correlation 
based  technique  [1].  The  output  of  this  pathway,  Tj(t),  is  taken 
to  be  an  estimate  of  the  intrinsic  torque. 

An  iterative,  cross-correlation  based  technique  [2]  is  then 
used  to  fit  a  Hammerstein  cascade  between  9(t)  and  the  resid¬ 
uals,  Tq(t)  —  where  77(f),  is  the  output  of  the  estimated 

intrinsic  dynamics. 

Finally,  note  that  the  reflex  torque  will  act  as  noise  in  the 
estimation  of  the  intrinsic  dynamics.  Thus,  it  may  be  possible 
to  improve  the  estimate  of  the  intrinsic  dynamics,  by  fitting  a 
new  IRF  between  9(t )  and  the  residuals  Tq(t)  —  TR(t).  This 
iteration  continues,  alternating  between  the  intrinsic  and  reflex 
dynamics,  until  the  model  accuracy  converges. 


B.  Separable  Least  Squares  Identification 

Parametric  optimization  methods  can  also  be  used  to  identi¬ 
fy  the  elements  of  the  PCSM.  Thus,  the  objective  would  be  to 
find  the  parameter  vector, 

P  =  [  hj3  h^s  cT  ]  (5) 

that  minimizes  the  cost  function 

l  N  2 

Vn(J3)  =  ^J2  W  -  Tq(t,  P))  (6) 

v  t=i 


where  Tq(t,  (3)  is  the  model  output  computed  using  (4),  where 
the  values  for  the  IRFs  and  polynomial  coefficients  are  con¬ 
tained  in  the  parameter  vector,  /?,  defined  in  (5). 

In  principle,  one  could  use  a  gradient  descent  procedure, 
such  as  the  Levenberg-Marquardt  (L-M)  algorithm,  to  solve 
the  minimization.  Thus,  one  would  start  with  an  initial  esti¬ 
mate,  and  refine  it,  in  the  case  of  L-M,  as  follows, 

/?(*+!)  =  ft®  +  (JTJ  +  Sk!)-1  JTe  (7) 


where  5k  is  a  regularization  parameter  used  to  control  the  con¬ 
vergence  rate  and  stability.  The  matrix  J  is  the  Jacobian,  a 
matrix  whose  [f ,  m]  entry  contains  the  partial  derivative  of  the 
model  output  at  time  t  with  respect  to  the  m’th  parameter. 
Thus, 


J(f,  to) 


dfq(t,p) 

dp(m) 


(8) 


However,  the  PCSM  has  (2 Lj  +  1)  +  (LR  +  1)  +  (Q  +  1) 
parameters,  corresponding  to  the  weights  in  the  intrinsic  and 
reflex  IRFs  and  the  polynomial  coefficients  that  describe  the 
nonlinearity.  Depending  on  the  memory  length  of  the  IRFs, 
and  the  order  of  the  polynomial,  the  stiffness  model  could  eas¬ 
ily  have  more  than  100  parameters!  Thus,  it  appears  that  the 
parametric  optimization  will  have  to  search  over  a  100+  di¬ 
mensional  parameter  space.  Furthermore,  each  step  (7)  in  the 
optimization  would  require  the  generation  and  inversion  of  a 
100  by  100  (or  larger)  matrix.  This  is  clearly  impractical. 

Note,  however,  that  if  the  nonlinearity  is  fixed,  the  output  of 
the  stiffness  model  (4)  is  a  linear  function  of  the  filter  weights 
of  the  two  IRFs.  Thus,  for  a  given  static  nonlinearity,  the  op¬ 
timal  weights  for  both  IRFs  can  be  found  simultaneously  in 
closed  form  by  solving  a  linear  regression. 

Subdivide  the  parameter  vector  into  linear  and  nonlinear  pa¬ 
rameters. 


P=[PT\Pl]=[hT 


IS 


hT 


RS 


(9) 


Then,  for  any  given  nonlinear  parameters  (polynomial  coeffi¬ 
cients),  pn,  the  optimal  linear  parameters  (IRF  weights),  Pi, 
can  be  found  by  solving  the  linear  regression. 


Tq  =  X(pn)pl  +  e  (10) 

where  Tq  is  a  vector  whose  t’th  element  is  Tq(t),  e  is  a  vector 
containing  the  errors,  and  the  regression  matrix,  X(pn),  con¬ 
tains  both  advanced  and  lagged  copies  of  9  (since  the  intrinsic 
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dynamics  are  two-sided),  and  lagged  copies  of  the  nonlinearity 
output.  Clearly,  only  the  latter  columns  depend  on  the  choice 
of  polynomial  coefficients. 

By  solving  (10),  (3i  is  expressed  as  a  function  of  the  nonlin¬ 
ear  parameters.  Thus,  the  model  output  and  cost  function  can 
be  considered  to  be  functions  of  the  nonlinear  parameters  only, 
and  written  as  Tq(t,  (3n)  and  Vat(/3„).  As  a  result,  the  gradient 
based  optimization  need  only  search  for  the  optimal  (3n,  reduc¬ 
ing  the  dimension  of  the  search  space  to  Q  + 1.  This  procedure 
is  known  as  separable  least  squares  optimization  [8], 

To  apply  L-M  to  the  SLS  problem,  the  separated  Jacobian 
must  be  computed.  First  consider  the  Jacobian  of  the  unsepa¬ 
rated  problem,  J,  which  can  be  computed  using  the  chain  rule. 
Partition  it  into  linear  and  nonlinear  columns,  J  =  [  Ji  |  Jn  ] 
corresponding  to  the  linear  and  nonlinear  parameters.  Then,  let 
Pi  be  an  orthogonal  projection  onto  the  column-space  of  J/.  It 
can  be  shown  [8]  that  the  Jacobian  for  the  separated  problem 
is  given  by, 

Jsls  =  (/  -  P,)Jn  (11) 

Thus,  to  use  the  L-M  optimization  method  on  a  SLS  problem, 
the  nonlinear  parameter  vector,  /3„,  is  updated  according  to  (7), 
but  using  Jsis  from  (1 1)  as  the  Jacobian. 


Typical  Experimental  Data 


Fig.  2.  Extract  from  a  typical  experimental  trial  showing  4  seconds  of  position, 
computed  velocity,  torque  and  GS  EMG. 


III.  Experimental  Results 

To  evaluate  the  utility  of  the  SLS  algorithm,  we  used  it 
to  estimate  intrinsic  and  reflex  stiffness  dynamics  in  a  spas¬ 
tic  spinal  cord  injured  (SCI)  patient.  Previously,  an  iterative, 
cross-correlation  based  algorithm  has  been  used  to  fit  the  PC- 
SM  to  experimental  data  [5],  [6].  Tests  based  on  this  model  are 
being  investigated  for  possible  clinical  use  [7].  Consequently, 
methods  for  efficient,  unbiased  estimates  of  its  elements  are 
potentially  very  important. 

The  experimental  paradigm  has  been  described  elsewhere 
[5].  Briefly,  the  subjects  had  an  incomplete  loss  of  motor  func¬ 
tion,  and  clinically  evident  spasticity  due  to  a  previous  spinal 
cord  injury[7],  which  was  associated  with  hyper-active  stretch 
reflexes.  They  lay  supine  with  their  left  foot  attached  to  a  ro¬ 


tary  hydraulic  actuator  by  a  custom  fitted  fiber-glass  boot.  An¬ 
kle  torque  and  position,  and  the  EMG  over  the  Gastrocnemius- 
Soleus  muscles  were  recorded.  The  subject  maintained  a  con¬ 
stant  background  contraction,  while  a  broad-band  position  per¬ 
turbation,  whose  spectrum  was  shaped  to  preserve  the  stretch 
reflex  [5],  was  applied.  Fig.  2  shows  4  seconds  of  the  30  sec¬ 
ond  (6000  points  at  200  Hz)  data  records  used  in  the  analysis. 


Parallel  Cascade  Stiffness  Model:  Correlation  Method 
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Fig.  3.  Elements  of  the  stiffness  model  identified  using  the  traditional,  iterative 
correlation  based  algorithm. 


We  examined  the  dynamic  relationship  between  the  ankle 
position  and  torque,  by  fitting  a  PCSM,  as  shown  in  Fig.  1, 
between  the  two  signals.  Models  were  fitted  between  the  first 
5000  points  of  the  measured  position  and  torque  signals.  The 
accuracy  of  the  identified  models  was  then  tested  by  using 
them  to  predict  the  ankle  torque  measured  in  the  remaining 
1000  points,  using  only  the  measured  position.  In  the  sequel, 
the  initial  5000  points  of  data  will  be  referred  to  as  the  identi¬ 
fication  segment,  whereas  the  final  1000  points  will  be  called 
the  validation  segment. 

First,  we  used  the  conventional,  iterative,  correlation  based 
algorithm  [5]  to  fit  the  intrinsic  and  reflex  pathways  of  the 
model.  Typical  values  were  chosen  for  the  IRF  lengths  and 
polynomial  order.  Thus,  we  set  the  length  of  the  intrinsic  dy¬ 
namics  to  be  ±35  ms,  the  length  of  the  reflex  dynamics  to  be 
320  ms,  and  the  polynomial  order  to  be  6.  Improvement  ceased 
after  3  iterations,  producing  the  model  shown  in  Fig.  3,  which 
fit  the  identification  segment  with  83.4  %  variance  accounted 
for  (%  VAF),  and  predicted  the  torque  in  the  validation  seg¬ 
ment  with  9 1 .2  %  VAF. 

Next,  we  used  the  separable  least  squares  approach  to  identi¬ 
fy  the  same  model  structure.  The  nonlinearity  identified  by  the 
iterative  correlation  based  algorithm  was  used  to  initialize  the 
SLS  procedure.  The  resulting  model  is  shown  in  Fig.  4.  It  has 
exactly  the  same  structure  as  the  model  identified  by  the  itera¬ 
tive  correlation  based  scheme,  but  accounted  for  86  %  VAF  in 
the  identification  segment,  and  93.7  %  VAF  in  the  validation 
segment,  both  significant  improvements  over  the  predictions 
made  by  the  first  model. 

In  comparing  the  two  identified  models,  the  IRFs  identified 
by  the  conventional  scheme  include  high  frequency  oscillatory 
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Fig.  5.  Errors  in  the  torque  predictions  during  2  sec.  of  the  validation  segment. 


Fig.  4.  Elements  of  the  stiffness  model  identified  using  the  separable  least 

squares  optimization  based  algorithm. 

components.  While  these  are  still  evident  in  the  SLS  mod¬ 
el,  they  are  greatly  reduced  in  amplitude.  The  nonlinearity 
identified  by  the  SLS  scheme  resembles  the  combination  of 
a  small  positive  threshold,  at  about  0.8  rad/sec,  followed  by  a 
half-wave  rectifier,  an  interpretation  consistent  with  empirical 
studies  of  the  stretch  reflex  [3].  The  functional  significance 
of  the  nonlinearity  identified  by  the  conventional  identification 
scheme  is  not  so  evident. 

Figure  5  shows  the  prediction  errors  due  to  the  two  identified 
models  for  2  seconds  of  the  validation  segment.  The  solid  line 
shows  the  residuals  due  to  the  model  identified  using  the  SLS 
technique,  whereas  the  dash-dotted  line  shows  the  errors  in  the 
output  of  the  model  identified  using  the  traditional  approach. 
The  improvement  in  prediction  accuracy  due  to  the  SLS  model 
is  clearly  evident  in  this  figure. 

IV.  Discussion 

In  this  paper,  we  presented  a  block  structured  model  for  the 
dynamic  stiffness  of  the  human  ankle,  and  developed  a  SLS 
identification  technique  for  it.  This  algorithm  was  tested  on 
experimental  data  from  a  SCI  patient,  and  the  resulting  mod¬ 
el  compared  to  one  identified  using  currently  available  tech¬ 
niques.  The  SLS  model  produced  better  predictions  of  the 
measured  torque,  both  in  the  training  sample,  and  on  validation 
data.  Furthermore,  the  elements  of  the  SLS  model  were  visibly 
less  noisy,  facilitating  their  physiological  interpretation. 

This  SLS  approach  is  well  suited  to  the  parallel  cascade 
stiffness  model  structure,  because  the  model  output  is  linear 
in  most  of  its  parameters.  Thus,  using  the  SLS  approach,  it  is 
possible  to  reduce  the  dimension  of  the  search  space  dramati¬ 
cally.  In  this  case,  the  model  had  102  parameters  (31  weights 
for  his(r),  64  weights  for  hRs(r),  and  7  polynomial  coeffi¬ 
cients  for  nns{'))-  However,  by  using  the  SLS  algorithm,  it 
was  only  necessary  to  search  over  a  7  dimensional  parameter 
space,  corresponding  to  the  7  polynomial  coefficients. 

The  biggest  disadvantage  with  using  a  gradient  descent 
based  algorithm,  is  that  the  iteration  may  converge  to  a  sub- 
optimal  local  minimum.  When  using  SLS,  it  is  essential  to 


find  an  initial  guess  at  the  nonlinear  parameters  which  is  close 
enough  to  the  globally  optimal  solution  that  the  search  will 
find  it,  instead  of  converging  to  a  sub-optimal  solution.  In  this 
study,  we  used  the  results  obtained  by  the  conventional,  cross¬ 
correlation  based  method  to  provide  this  initial  guess. 

The  static  nonlinearity  identified  by  the  SLS  method  includ¬ 
ed  a  small  positive  deflection  at  relatively  large  negative  veloc¬ 
ities.  Given  the  underlying  physiology,  this  is  likely  to  be  an 
artifact,  probably  due  to  the  use  of  polynomials  in  represent¬ 
ing  the  nonlinearity.  Other  representations,  including  rational 
polynomials,  splines,  or  sigmoidal  neural  networks,  may  be 
more  appropriate.  The  incorporation  of  these  nonlinearities  in¬ 
to  the  SLS  algorithm  is  a  topic  of  ongoing  research. 
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